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Intergranular fracture is a dominant mode of failure in ultrafine grained materials. In the present 
study, the atomistic mechanisms of grain-boundary debonding during intergranular fracture in 
aluminum are modeled using a coupled molecular dynamics - finite element simulation. Using a 
statistical mechanics approach, a cohesive-zone law in the form of a traction-displacement 
constitutive relationship, characterizing the load transfer across the plane of a growing edge 
crack, is extracted from atomistic simulations and then recast in a form suitable for inclusion 
within a continuum finite element model. The cohesive-zone law derived by the presented 
technique is free of finite size effects and is statistically representative for describing the 
interfacial debonding of a grain boundary (GB) interface examined at atomic length scales. By 
incorporating the cohesive-zone law in cohesive-zone finite elements, the debonding of a GB 
interface can be simulated in a coupled continuum-atomistic model, in which a crack starts in the 
continuum environment, smoothly penetrates the continuum- atomistic interface, and continues 
its propagation in the atomistic environment. This study is a step towards relating atomistically 
derived decohesion laws to macroscopic predictions of fracture and constructing multiscale 
models for nanocrystalline and ultrafine grained materials. 
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1. Introduction 

Physics-based modeling of fracture begins at nanometer dimensional scales in which atomistic 
simulation is used to predict the evolution of various deformation mechanisms on a fundamental 
level. These mechanisms include dislocation nucleation and interaction, interstitial void 
formation, and atomic diffusion. The development of these damage mechanisms progresses into 
microscale processes such as local plasticity and small crack formation. Ultimately, damage 
progression leads to macroscopic failure modes such as plastic yielding and large cracks 
exhibiting Mode I, II, and III behavior. To capture this macroscopic behavior, the use of direct 
atomistic simulations is impractical due to their strong inherent limitations in size and time 
scales. Thus, multiscale modeling strategies are needed for the purpose of minimizing 
computational requirements and developing a unified description of the hierarchy of material 
processes that govern fracture. 

Recently, a multiscale strategy using a constitutive relation-based continuum modeling of the 
process of intergranular fracture in aluminum has been proposed [1,2]. In this strategy, 
information from the lower atomistic level is transferred to the upper continuum level through 
constitutive relations governing different damage modes. The constitutive relations must reflect 
the most dominant characteristics of the individual deformation mechanisms revealed at the 
lower atomistic level to adequately model the behavior observed at higher length scales. For 
fracture simulations, the traction-displacement relationship (as introduced by Dugdale and 
Barenblatt [3,4]), is used to model the cohesive properties of the grain boundary (GB) interface 
for a specific fracture mode. The traction-displacement relationship in this strategy is extracted 
from atomistic molecular dynamics (MD) simulations [5], and is incorporated into cohesive-zone 
models (CZM) [6] used in finite element models [7]. This yields a finite element analysis 
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capability to study the behavior of metallic microstructural damage at the grain scale [1,2]. In 
this multiscale approach, the CZM elements represent the mechanical response of the GBs as 
obtained from atomistic simulations, while the constitutive properties of the crystalline grains are 
represented by continuum finite elements with appropriate anisotropic elastic properties. Thus, 
an integrated multiscale modeling strategy for understanding the connection between 
macroscopic fracture and the underlying atomistic failure mechanisms emerges. In principle, the 
range of relevant length scales achieved by this strategy can span from the atomistic (10" 9 m) to 
the macroscopic (>10‘ 3 m). 

An essential part of this strategy is the extraction of CZMs for interface decohesion from an 
atomistic simulation. Efforts on this topic have been made by various groups in the last few years 
[8-11]. All of those simulations showed a highly elevated stress of debonding ranging from 15 to 
50 GPa, which are generally two orders of magnitude higher than the experimentally observed 
strength of the corresponding material [8], There are various reasons for this discrepancy. For 
example, the approach used in all of these works [8-11] was based on simulating the debonding 
of a flat interface under a constant tensile strain rate perpendicular to the interface. The dynamics 
of the atoms was severely constrained by the boundary conditions, which did not allow for 
Poisson contraction and prevented shear deformation. As a result, plastic processes, such as 
dislocation slip, interface sliding and interface diffusion, were strongly suppressed. 
Consequently, the simulated mechanism for interface decohesion in these works reproduced an 
idealized process of atomic adhesion (strength) rather than fracture at the interface. 

Recently, a methodology for extracting CZMs from atomistic simulations of crack 
propagation has been developed [5]. The main goal of this methodology was to extract and 
understand the contributions of different atomistic processes to an MD-based CZM decohesion 
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law for intergranular fracture under local conditions of a propagating crack. A new concept of 
defining cohesive-zone-volume elements (CZVE) as atomic subdomains encompassing the 
debonding interface, and presenting an atomistic equivalent of CZM elements was introduced. In 
this concept, the CZM is viewed as a statistical representation of a large, Gibbs-type ensemble of 
CZVEs placed along the crack path during crack propagation. The CZM is obtained as a 
statistical average of the traction-displacement states of the CZVE ensemble under the conditions 
of a steady-state crack propagation. 

In this paper, the concept of implementing a statistical CZVE ensemble to extract a CZM 
constitutive law in terms of a traction-displacement relationship will be used in an improved 
coupled atomistic-continuum model for intergranular fracture in aluminum. The dependence of 
the loading conditions and temperature on the CZM will be determined. The derived CZM will 
be incorporated within finite elements introduced in the continuum domain of the coupled 
atomistic-continuum model as an extension of the atomistic interface. This provides the ability to 
simulate crack growth smoothly propagating from the continuum region into the atomistic 
domain. This technique allows multiscale atomistic-continuum modeling of fracture processes 
where fracture may not be confined to the atomistic domain. 

2. The simulation approach 

The simulation approach used in this study is based on a coupled atomistic-continuum model 

that allows a large atomistic domain (containing 10 2 * * * 6 or more atoms) to be embedded within a 
continuum domain of micron dimensions. The atomistic domain is simulated using MD, where 

the interatomic forces are represented by an atomistic potential suitably fitted to reproduce the 

material properties of aluminum [12]. The continuum domain is simulated by using the finite 
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element method (FEM) with anisotropic elastic properties derived from the aluminum potential 
[12] at room temperature (300K). The coupling between FEM and MD is achieved by the 
recently-developed embedded statistical coupling method (ESCM) [13,14] to provide elastic 
boundary conditions for the atomistic domain and transfer the applied far field mechanical load 
into the atomistic system. 

The ESCM method is based on the concept that the continuum representation of a material, 
in terms of stress-strain fields, is a statistical representation of its atomic structure. Following this 
idea, the connection between the atomistic and the continuum representations at the MD/FE 
interface is performed using a statistical mechanics approach. The method uses statistical 
averaging over both time and volume of atomistic subdomains at the MD/FE interface to provide 
nodal displacement boundary conditions to the continuum FEM domain. The FEM generates 
interface reaction forces that are uniformly distributed over the interface atoms in the form of 
constant traction boundary conditions [15] to the MD domain. Thus, this approach is based on 
solving the special boundary value problem (BVP) at the MD/FE interface for a coupled MD- 
FEM system and may be described as a local-nonlocal BVP because it relates local continuum 
nodal quantities with statistical averages of nonlocal atomistic quantities over selected atomic 
volumes. For this study, one finite element at the interface encompasses a region of several 
hundred to several thousand atoms, which positions have been averaged over 1 ps time interval. 
At those scales, the discreteness of the atomic structure is sufficiently homogenized so that the 
FEM domain responds to the atomistic domain as an extension of the continuum. The constant 
traction boundary conditions of the MD domain ensure that the elastic field from the FEM 
domain is correctly transferred to the atomistic region. An iterative procedure using the MD 
statistical displacements establishes a balance between the FEM-computed forces and the MD- 
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computed forces at the interface. This force balance ensures stress-displacement continuity at the 
interface. 

The model used to extract CZMs for a given grain boundary in aluminum is presented in Fig. 
1 . A circular MD domain with a diameter duo = 90 nm was embedded in a square FEM mesh 
with side dimension c/fe = 20 djyiD- In this example, the system represents a bicrystal with a <1 1 
0> X99 symmetric tilt GB passing along the line (y = 0) in Fig. 1. In the MD domain, the grain 
boundary was formed by symmetrically rotating the two crystals at 89.42° with respect to each 
other around a common tilt axis coinciding with the [110] crystallographic direction, chosen as 
the z-direction in Fig. 1. The system thickness in the z- direction was h = 2.9 nm and includes 
10(1 1 0) crystallographic stacking planes. Periodic boundary conditions were applied in the z- 
direction of the MD system to emulate a bulk atomic state in this direction. Structurally, the 
atomistic model is equivalent to the one discussed in detail in [5]. The MD stress along the 
periodic z-direction in the MD domain was maintained at zero by applying the Parrinello- 
Rahman constant stress technique [16]. The MD simulation was performed at constant 
temperature using the Nose-Hoover thermostat [17]. 

3. Atomistic derivation of CZM 

To derive the CZM for interface debonding, an edge crack along the GB line (Fig. 1), which 
extended from one of the outer boundaries of the FEM mesh to the interface with the MD 
domain, was inserted in the coupled MD-FEM model. The crack was propagated into the MD 
domain by imposing a tensile uniaxial strain of 1.5% in the y-direction (perpendicular to the 
crack plane). The atomistic processes associated with crack propagation along this GB were 
studied in detail in [5], Energetic analysis of the different dissipative mechanisms during crack 
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growth was given in [18], and some aspects of the dynamics of the fracture process were 
reported in [19]. 

CZVEs of size / x 2/ x h in the x-, y- and z- directions, respectively, with / = 1.9 nm, were 
introduced along the GB plane as discussed in [5] (Fig. 2). Each CZVE contained approximately 
1250 atoms over which the stress-displacement response of the element was calculated. During 
crack propagation, the data from each CZVE was used to extract profiles of the normal stress, 
cr (x), and opening-displacement, X(x) , along the GB interface. Figure 2 shows an example of 
one of these profiles, together with a snapshot of the atomistic configuration, from which it was 
extracted. A fit, using o yy = K,/ ^2jt(x-x 0 ) , gave an estimate of the stress intensity factor of K\ 
= 0.88 [MPa-m ] for x 0 defined at the position of the crack tip. This value is slightly lower than 
the stress-intensity factor for a pure linear elastic model of the same geometry, K\ = 1.09 [MPa- 
m ■], estimated from [20]. The discrepancy is most likely a result of the developed plasticity at 
the crack tip in the atomistic simulation. 

Using a statistical mechanics approach [21], the CZM is derived in the following way. The 
state of each CZVE at position x and instant of time t is defined by the state variables a (x,i) 

and X(x,t). The CZVE state is represented as a point ( X, a vv ) in a cr yy -vs-X configuration space. 
In the statistical limit of steady-state crack propagation that occurs over an infinitely long time 
over an infinitely long interface, all realized CZVE states would produce a p(X,a vy ) distribution 

(density of states) that is a continuous function independent of time. The CZM is defined as a 
statistical average, a (X), of a over a small interval (X - AX, X + AX) of this distribution 
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X+AX 00 

oM = (°M) = (i) 

f f p( x '’°v,K a ' 

X'=X-AX a^,=- oo 


at the limit of AX -» 0 . 

Numerically, p(X,a yy ) can be approximated by using the already described coupled MD- 
FEM simulation model. For a finite number of measurements, p(X,cr yy ) takes a discrete fonn 

P( K ' a yy ) = 2 8 ( X - K / ) 8 K “ a xv ( X / )) ’ ( 2 ) 

i 


where the sum is over all measured data points (X ( ,a vv (X ( )j and 6(.r-x ; ) is the Dirac delta 
function. 

After substituting p(X,a vv ) from Eq. (2) into Eq. (1), a vv (>.) takes the form of a moving 
average 




N^Kj E |X — AX, X + AXjj ■ 




( 3 ) 


where the sum is taken over all N points (X,a vv ), which are between X- AX and X + AX. The 

choice of the size of the interval, AX , depends on the desired resolution of the CZM with respect 
to X. 

Figure 3 shows the derived p(X,a yy ) for the example given in Fig. 2 (the case for K\ {) = 1.09 

[MPa-m "]), after collecting data from 180 snapshots that were taken at intervals of 1 ps. 
Another example with a lower stress intensity of K\ 2) = 0.72 [MPa-m 1 2 ] with all other conditions 
being equal, is also given for comparison. The solid lines, drawn through both distributions in 
Fig. 3, are the averaged a (X) CZM curves (Eq. 3). Both CZMs were derived using a resolution 


of AX = 0.1 mn. 
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The comparison between the two CZMs in Fig. 3 indicates that increasing the stress intensity 
in the simulation increases the displacement at full debonding, Ao (1) > Ao <2) , while the peak stress, 
a p = 6.4 GPa, is approximately constant. Here, the increase in K\, is associated with an increase 
in plastic mechanisms (i.e., increased nucleation of dislocations) in the softening region rather 
than an increase in the force needed to break atomic bonds. 

A substantial effect of temperature can be seen on the CZM result at T=100K and K\ = 0.72 
[MPa-m ‘] shown in Fig. 3 when compared to the CZM result at T=300K shown in Fig. 4 for the 
same loading conditions. Increasing the temperature reduced the peak stress from 6.3 GPa to 5.4 
GPa, but increased the displacement at full debonding from 1.7 nm to 2.7 mn. The implication is 
that the work, T = J a vr (/,)d'k (Fig. 4), done in the softening region increases due to increased 

crack tip plasticity. The decrease in a p can be explained by the decrease of force required for 
bond breaking due to the increased atomic thermal motion at higher temperature. 

4. Crack propagation from a continuum to atomistic region 

The next logical step after extracting the CZM curve from atomistic simulations is to 
implement it within a CZM element and verify that the debonding predicted by the atomistic 
simulation is reproduced by the finite element analysis. For this purpose, the same MD-FEM 
coupled model as described in Section 2 was used with a slight, but essential, modification: the 
atomistic GB in the MD domain was extended into the continuum domain by placing CZM 
elements along the interface between the two grains. The thickness of the system was the same 
as that used for extracting CZMs, h = 2.9 nm, while the lateral dimensions were reduced by a 
factor of two to save simulation time such that duo = 45 nm and d\-\- = 900 nm. It was verified 
that decreasing of the system size did not affect the extracted CZM. 
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Again, an edge crack was inserted along the GB line, but unlike the model in Section 3, the 
crack did not extend all the way to the MD domain. Rather, it had to propagate through four 
CZM elements (12 nm) in the continuum before reaching the MD/FE interface. Its growth into 
the continuum domain is governed by the cohesive zone model incorporated within the CZM. 
For computational efficiency, this CZM law was defined as a least-square bilinear fit to the 
extracted p(X,a vv ) distribution at T = 300 K, shown in Fig. 4. The two points, (k ;) , a ;) j and 
(X 0 ,0) were used as fitting parameters (Fig. 5). The so-called “penalty stiffness”, defined by the 
slope of the CZM curve before the maximum traction is reached, was set to match the initial 
slope (from a = 0 to a = a p ) of the atomistic traction-displacement relationship. In this way, the 
stiffness of the CZM before opening matched the elastic stiffness of the GB interface. 

A uniform tensile strain that gradually increased from 0 to 1.5% was applied to the outer 
boundaries of the FEM domain. The strain was applied in 0.1% increments over 10 ps MD time 
intervals (equivalent to 1x10 s s' 1 strain rate). This constant strain rate allowed a gradual growth 
of the edge crack, first in the continuum domain, and then penetrating into the MD domain. The 
process is shown in several sequential snapshots in Figs. 6(a-f) that focus on the MD domain in 
the center of the coupled MD-FEM model. In Figs 6(a) and 6(b), the crack propagates through 
the FEM mesh, sequentially opening the CZM elements placed along its path. In Fig. 6(c), the 
crack tip enters the MD domain. The snapshot depicts the decohesion at the MD/FE interface in 
unison with the opening of the last CZM element. In Figs. 6(d) and 6(e), the crack continues its 
growth in the MD domain. Plastic processes start to develop in the form of dislocation 
nucleations, seen as short green lines marking the stacking faults at the extended dislocation 
cores on both sides of the crack tip. Fig. 6(f) presents a crack that has propagated well inside the 
MD domain. Active dislocation nucleation accompanies the propagation. The snapshots in Figs 
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6(a-f) also quantitatively show that both the crack opening and the stress-displacement field are 
continuous at the MD/FE interface during the entire process of crack propagation. 

5. Summary 

This study examined a methodology for the atomistic derivation of a cohesive-zone law that 
can be implemented in cohesive-zone finite element models for simulating fracture in 
nanocrystalline or ultrafine grained aluminum. The cohesive-zone law is derived as a statistical 
mechanics representation of crack propagation through defining an ensemble of atomistic 
cohesive-zone-volume-elements. The simulation methodology implements a coupled MD-FEM 
model utilizing statistical coupling between the atomistic and the continuum representations of 
the material. The model allows large atomistic domains (containing 10 6 or more atoms) 
simulated by MD at finite temperature to be embedded into an FEM mesh of micrometer 
dimensions. The cohesive-zone laws derived by this technique are free of finite size effects and 
are statistically representative for describing the interfacial debonding of an idealized interface 
examined at atomic length scales. As a test example, a coupled continuum-atomistic model was 
presented, in which a crack starts in the continuum environment, smoothly penetrates the 
continuum-atomistic interface, and continues its propagation in the atomistic environment. The 
model introduces the possibility of conducting multiscale simulations where fracture is not 
strictly limited to the atomistic region, which may greatly benefit the development of multiscale 
atomistic-continuum models of nanocrystalline and ultrafine grained materials. 
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Fig. 1 . Schematic description of the ESCM coupled molecular-dynamics - finite-element model 
configuration. 
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Fig. 2. Stress and opening profiles from cohesive-zone-volume-elements (CZVE) along a grain- 
boundary edge crack growing inside the molecular-dynamics region. 
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Fig. 3. Statistical distributions of measurements from individual CZVEs (dots), and averaged 
traction-displacement relationships (solid lines) at T = 100 K extracted from the ESCM model at 
different stress intensity factors. The values a p , A(/ ' \ and Xo <2) , indicate the peak stress, and the 
displacements at full debonding at Ki (1) and Ki <2) , respectively. 


17 



Contributed paper proceeding to the 2008 TMS meeting, New Orleans, LA, March 9-13, 2008. 




V 10 °) 


8 

6 

4 

2 

0 


y 30 °) 




0 1 

X [nm] 


a ( 10 °) 
p 


a 


T = 1 00 K 
— T = 300 K 


Fig. 4. Traction-displacement relationships o yy (k), and energy of decohesion T(X), (the solid 
lines) at T = 100 and 300K. The statistical distributions of measurements from individual 
CZVEs, used to extract the traction-displacement relationships are also given (dots) for 
comparison. The values a p , <100) , X 0 <100) , a p <300) , and X 0 <300) , indicate the peak stress and the 
displacement at full debonding at T = 100K and 300K, respectively. 
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Fig. 5. Bilinear CZM, parameterized as a least-square fit to the extracted statistical stress- 
displacement distribution of measurements from individual CZVEs at T = 300K. The values o p , 
X p , and Xo are used as parameters. 
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Fig. 6. Structural snapshots with imposed a yy stress map on the f.c.c atoms of the coupled MD- 
FEM/CZM model simulating an edge crack passing through the MD-FEM interface and 
debonding a GB in aluminum. Strain rate of 1x10 s s’ 1 was applied. Atoms in h.c.p state 
(indicating stacking faults produced from partial dislocations) are colored in green, and the rest 
(disordered and surface atoms) are shown in black. The six snapshots are taken at uniaxial tensile 
strain £yy as follows: (a) 1.2%; (b) 1.3%; (c) 1.3%; (d) 1.4%; (e) 1.5%; (f) 1.5%. 
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